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A  mathematical  model  of  a  single  unit  solid  oxide  fuel  cell  (SOFC)  based  on  partial  differential  equations 
is  presented  in  this  paper.  The  model  consists  of  mass,  momentum  and  charge  conservation  equations, 
relating  gas  concentration  to  voltage,  current  density,  and  other  relevant  fuel  cell  parameters.  The  model 
is  applied  and  validated  on  a  planar  disk-shape  anode-supported  cell.  Data  regarding  materials’  prop¬ 
erties  and  kinetic  parameters  are  taken  from  the  open  literature  or  derived  from  experimental  results. 
A  commercial  tool  using  the  finite  element  method  is  used  to  solve  the  set  of  equations.  The  simulated 
steady  state  performance  of  the  fuel  cell,  for  different  operating  conditions  and  current  densities,  is  further 
investigated. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  structure  of  a  mathematical  model  for  fuel  cells  (FC)  simu¬ 
lation  usually  depends  on  the  purpose  of  the  model  itself.  When, 
for  example,  the  main  purpose  of  the  model  is  to  provide  FC  per¬ 
formance,  in  order  to  analyze  the  whole  system  (or  power  plant) 
performance,  the  variation  of  physical-chemical  variables  (such  as 
gas  concentration,  temperature,  pressure  and  current  density,  for 
example)  are  not  relevant;  however  the  performance,  in  terms  of 
power,  heat  and  input  requirements  are  important.  On  the  other 
hand,  if  the  focus  of  the  investigation  is  the  design  of  a  stack,  a 
single  cell,  or  even  a  specific  component  of  the  cell,  a  model  that 
considers  a  detailed  description  of  each  phenomenon  taking  place 
within  the  fuel  cell  is  needed. 

According  to  the  level  of  details  considered,  mathematical  mod¬ 
els  for  fuel  cells  may  be  classified  as: 

•  Black-box  models  [  1  ],  where  single  cells  or  stacks  are  represented 
as  non-dimensional  devices  with  a  number  of  mass  and  energy 
inputs  and  outputs.  Transport  phenomena  occurring  within  the 
system  are  not  simulated  with  this  approach,  but  mass  and  energy 
balances  are  solved  for  computing  the  inlet  and  outlet  streams’ 
properties,  as  well  as  for  estimating  fuel  cell  power  and  efficiency. 
These  models  are  not  suitable  to  analyze  phenomena  occurring 
within  a  fuel  cell,  but  they  are  a  necessary  tool  for  analyzing  how 
the  fuel  cell  affects  other  parts  of  the  system  [2,3]. 
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•  Single  cell  models  with  details  of  the  cell  components  [4-10]. 
These  models  describe  the  behaviour  of  concentration  of 
chemical  species,  gas  velocity,  temperature,  current  and  the  dis¬ 
tribution  of  other  physical  parameters  in  the  fuel  cell,  and  they 
are  usually  used  as  design  tools  to  optimize  geometry  parameters 
and  dimensions  of  the  single  cell  or  of  the  stack  (e.g.  electrodes 
thickness,  rib  width)  and  to  predict  the  operating  conditions. 

•  The  last  approach  is  the  microscopic  one,  in  which  diffusion  and 
adsorption  phenomena  are  studied  in  detail  and  parameters  used 
in  the  macroscopic  models,  such  as  porosity,  tortuosity,  and  acti¬ 
vation  energy  can  be  estimated. 

The  last  two  model  approaches  are  more  complex  than  the  first 
one,  and  numerical  methods  must  be  used  for  solving  the  systems 
of  partial  differential  equations. 

Mathematical  models  describing  electrochemical  and/or  fluid 
dynamic  performance  of  the  fuel  cell  may  be:  one-dimensional 
when  parameters  are  considered  to  vary  only  in  one  direction  (usu¬ 
ally  the  flow  direction)  [9],  two-dimensional  when  phenomena  are 
modelled  to  vary  over  a  cross-section  as  representative  behaviour 
of  the  entire  cell  performance  [4,5]  and  three-dimensional  when 
phenomena  are  considered  to  vary  within  the  whole  domain  and 
parameter  variations  in  every  direction  are  taken  into  account.  Sim¬ 
plifications  and  assumptions  which  lead  to  ID  and  2D  models  must 
be  carefully  chosen  to  minimize  the  distortion  of  reality  that  they 
introduce  and  to  avoid  neglecting  of  those  phenomena  which  are 
relevant  for  the  possible  application  of  the  model. 

Most  mathematical  models  available  in  the  open  literature 
do  not  consider  the  evolution  of  all  the  parameters  within  all 
the  dimensions  of  the  domain  considered.  Even  for  3D  mod¬ 
els,  the  most  frequent  approach  is  to  assume  a  lumped  Positive 
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Nomenclature 

Dij 

diffusion  coefficient  in  a  binary  system  [m2  s-1  ] 

Due 

Knudsen  diffusion  coefficient  [m2  s-1  ] 

A'.eff 

effective  diffusion  coefficient  [m2  s-1  ] 

e 

percentual  error  [-] 

Lact 

activation  energy  [J  mol-1  ] 

Eo 

ideal  Nernst  voltage  at  standard  pressure  [V] 

F 

Faraday  constant  [C  mol-1  ] 

j 

current  density  [Am-2] 

Jo 

exchange  current  density  [Am-2] 

^channel 

height  of  the  gas  channel  [m] 

Mi 

molecular  weight  of  component  i  [g  mol-1  ] 

M 

molecular  weight  of  gas  mixture  [g  mol-1  ] 

Ni 

mass  flow  of  component  i  [kgm-2  s-1  ] 

n 

number  of  electrons  involved  in  the  electrochemical 
reaction  [-] 

OCV 

open  circuit  voltage  [V] 

P 

pressure  [Pa] 

^ref 

reference  pressure  [Pa] 

Pi 

partial  pressure  of  component  i  [Pa] 

pi 

partial  pressure  of  component  i  at  the  elec¬ 
trode/electrolyte  interface  [Pa] 

Qan.cat 

volumetric  anode  or  cathode  flow  rate  [kg  s-1  ] 

tmanif 

radius  of  the  cell  section  where  the  flow  configura¬ 
tion  changes  from  normal  to  radial  [m] 

Rg 

universal  gas  constant  [J  mol-1 1<-1  ] 

T 

cell  temperature  [K] 

To 

standard  temperature  [K] 

u 

gas  velocity  [m  s-1  ] 

E^an.cat 

gas  inlet  velocity  [ms-1] 

Uf 

coefficient  of  utilization  [-] 

Vrev 

reversible  voltage  [V] 

'/cell 

cell  voltage  [V] 

Vi 

mass  fraction  of  component  i  [-] 

Subscripts 

act 

activation 

an 

anode  side 

cat 

cathode  side 

cone 

concentration 

el 

electric 

ion 

ionic 

Superscripts 

lab 

laboratory 

sim 

simulation 

Greek  letters 

&OCV 

corrective  coefficient  [-] 

P 

symmetry  factor  of  Butler-Volmer  equation  [-] 

£ 

porosity  [-] 

0 

potential  [V] 

y 

pre-exponential  factor  [Am-2] 

[i 

dynamic  viscosity  [kgm-1  s-1  ] 

T] 

voltage  loss  [V] 

P 

density  [kgm-3] 

Pel 

electronic  density  [Cm-3] 

P  ion 

ionic  density  [Cm-3] 

O 

electric/ionic  conductivity  [S  m-1  ] 

X 

tortuosity  [-] 

Electrode-Electrolyte-Negative  Electrode  (PEN)  structure  for  mod¬ 
elling  the  current  density  and  electrical  potential  distribution  in  the 
direction  perpendicular  to  the  flow  direction  [4-9]. 

The  result  of  such  assumption  is  a  fast  converging  solution.  On 
the  other  hand,  information  on  current  density,  electrical  potential 
distribution  and  gas  species  concentrations  within  the  PEN  domain 
is  neglected.  Furthermore,  such  models  are  restricted  to  a  specific 
geometry.  They  cannot  simulate  other  fuel  cells  with  slight  geom¬ 
etry  variations.  Bjaradwaj  et  al.  [4,5],  for  example,  developed  two 
different  models  for  simulating  the  operation  of  two  SOFCs,  one 
tubular  and  the  other  flat  tubular  (both  SOFC  types  are  produced 
by  Siemens  Westinghouse  Power  Corporation). 

The  time  consuming  process  of  defining  a  new  model  for  every 
geometry  modification  is  not  present  if  all  the  equations  are  in  a 
partial  differential  form  (PDE)  [11-16].  This  feature  is  of  partic¬ 
ular  importance  during  the  design  and  improvement  phase  of  a 
SOFC,  when  different  geometrical  configurations  need  to  be  com¬ 
pared. 

Bove  and  Ubertini  [16]  show  how  most  models  found  in  the  lit¬ 
erature  can  be  considered  simplifications  of  a  mathematical  model, 
entirely  based  on  partial  differential  equations,  where  the  variation 
of  all  the  involved  variables  is  considered  in  the  entire  domain. 
Limitations  and  simplifications  introduced  with  each  modelling 
approach  are  also  highlighted. 

In  the  present  study,  the  mathematical  model  described  in  [16] 
is  applied  and  validated  on  a  planar  disk-shape  anode-supported 
cell.  The  same  model  was  previously  applied  for  a  micro-tubular 
SOFC  geometry,  without  significant  modifications  [15].  The  results 
allowed  one  to  evaluate  different  performance  of  a  micro-tubular 
SOFC,  relative  to  different  current  collector  configurations. 

In  the  present  study,  the  model  is  presented  in  a  3D,  time- 
dependent  form,  and  is  implemented  for  a  2D  geometry,  under 
steady-state  conditions. 

2.  Fuel  cell  geometry 

The  model  is  applied  to  a  planar  disk-shape  SOFC  shown  in  Fig.  1 
currently  produced  by  InDEC  B.V,  and  operated  at  the  Fuel  Cell  Lab 
of  the  University  of  Perugia. 

The  anode  material  is  a  cermet  nickel  oxide  doped  with  yttrium 
stabilized  zirconia,  8%  yttrium  (NiO/8YSZ).  The  cathode  is  made 
of  8YSZ  with  strontium-doped  LaMn03  (8YSZ/LSM),  and  the  elec¬ 
trolyte  consists  of  dense  8YSZ. 

As  shown  in  Fig.  1 ,  the  anodic  and  cathodic  flows  are  introduced 
perpendicularly  to  the  plane  of  the  cell,  on  the  two  opposite  sides, 
through  two  alumina  tubes  (not  represented  in  Fig.  1).  The  anodic 
and  cathodic  flows  have  a  mirror  symmetric  configuration,  thus 
hereinafter  in  this  section,  the  description  of  the  anodic  or  cathodic 
flow  is  simply  referred  to  as  “flow”,  indifferently  if  referred  to  the 
anode  or  to  the  cathode. 

The  flow  reaches  the  electrode  surface  in  a  circular  region, 
where  it  changes  its  direction  from  perpendicular  to  radial,  with 
respect  to  the  fuel  cell  plane,  and  proceeds  along  the  electrode 
surface  in  a  co-flow  configuration.  The  electrical  current  is  drawn 
from  the  cell  through  a  mesh  (Pt  for  the  cathode  and  Ni  for  the 
anode),  located  between  the  cell  housing  and  the  electrode.  In 
the  present  study,  the  effect  of  such  mesh  on  the  velocity  field 
of  the  gas  in  the  gas  channel  is  neglected,  due  to  the  laminar 
regime  within  the  whole  cell.  Fig.  2  shows  the  Reynolds  num¬ 
ber  in  both  gas  channels  for  V=  0.656  V.  The  effect  of  the  different 
gas  flows  in  the  channels  on  Reynolds  number  values  is  clearly 
observed  since  cathodic  gas  flow  exceed  five  times  the  anodic 
one. 

The  geometrical  dimensions  of  the  single  cell  are  presented  in 
Table  1. 
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Oxidant  flow 
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Fig.  1.  Planar  disk-shape  SOFC  considered  in  the  present  study. 
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Fig.  2.  Reynolds  number  in  gas  channels  for  V=  0.656  V. 


3.  Process  description 

As  shown  in  Fig.  1,  once  the  gas  starts  flowing  radially,  it  is 
in  contact  with,  and  permeates,  the  porous  electrode  (anode  or 
cathode).  Diffusion  takes  place  in  the  whole  electrode  domain,  and 
proceeds  towards  the  reaction  zone  which  is  very  close  to  the  elec¬ 
trode/electrolyte  interface.  Since  both  electrodes  are  electrically 
conductive,  electrons  coming  from  (going  to)  an  external  circuit 
flow  through  it.  Oxygen  provided  at  the  cathode  side,  once  hav¬ 
ing  reached  the  reaction  zone,  combines  with  electrons  to  produce 
oxygen  ions,  according  to  the  following  reaction: 

io2+2e-^02-  (1) 

The  electrolyte  presents  high  electrical  resistance  and  high  ionic 
conductivity,  thus,  ideally,  only  ions  flow  through  it.  Moreover  the 


Table  1 

Tested  planar  disk-shape  dimensions. 


Active  area 

4848  mm2 

Anode  thickness 

545  |Jim 

Cathode  thickness 

40  pan 

Electrolyte  thickness 

5  fxm 

Cathode  channel  height 

1  mm 

Anode  channel  height 

0.6  mm 

Channels  width 

1.51  mm 

dense  structure  does  not  allow  any  molecular  permeation.  At  the 
anode  side,  ions  coming  from  the  electrolyte  react  with  hydrogen: 

H2  +  02“  ->■  H20  +  2e-  (2) 

Reactions  (1)  and  (2)  take  place  where  there  is  a  simultaneous 
presence  of  electrons,  ions  and  gases,  i.e.  in  the  so-called  triple¬ 
phase-boundary  (TPB).  As  a  consequence  of  the  reactions,  water, 
electricity  and  heat  are  produced.  Water  generated  in  the  anode 
side  diffuses  to  reach  the  gas  channel  and  leaves  the  cell.  Electrons 
released  in  reaction  (2),  migrating  through  the  anode,  are  collected 
through  the  current  collector,  and  are  transferred  to  the  cathode 
through  an  external  load.  At  the  cathode  side,  these  electrons  flow 
from  the  current  collector  to  the  cathodic  reaction  zone.  In  this  way, 
a  close  external  electrical  circuit  is  created. 

4.  Mathematical  model 

A  complete  mathematical  model  must  consider  all  the  transport 
phenomena  occurring  within  the  fuel  cell,  as  described  in  Section 
3.  Heat  transport,  charge  transport,  mass  diffusion  and  convection 
take  place  simultaneously  in  several  domains  of  the  fuel  cell. 

In  the  present  study,  each  sub-domain  (i.e.  anode,  cathode,  and 
electrolyte)  is  modelled  independently  from  the  other  following 
the  model  proposed  by  Bove  and  Ubertini  [16],  and  no  lumped 
structures  are  considered.  Phenomena  occurring  in  different 
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Fig.  3.  Velocity  field  in  a  longitudinal  section  of  the  anode  (15|jim  from 
anode/electrolyte  interface). 


10"4 


Fig.  4.  Velocity  field  in  a  longitudinal  section  of  the  cathode  (10|jim  from  cath¬ 
ode/electrolyte  interface). 


sub-domains  are  linked  to  each  other  by  appropriate  boundary  or 
interface  conditions. 

Energy  balance  is  not  included  because  the  SOFC  used  in  the 
experiments  was  operated  at  constant  temperature,  guaranteed  by 
an  electric  oven.  The  validation  of  the  energy  balance  would  thus 
be  extremely  hard  to  obtain,  and,  more  importantly,  the  isothermal 
conditions  would  not  reflect  the  real  operation  of  a  single  cell  of  a 
stack.  It  must  be  noted,  in  fact,  that  temperature  distribution  is  very 
non-uniform  within  a  single  cell  embedded  in  a  stack  [17],  while  in 
the  present  situation,  the  cell  is  mostly  at  isothermal  condition. 

Some  assumptions  have  been  made  to  simplify  conservation  and 
constitutive  laws: 

(1)  Electrochemical  reactions  are  confined  at  the  elec¬ 
trode/electrolyte  interface. 

(2)  Constant  temperature  within  the  cell 

(3)  Constitutive  law  for  ideal  gas 

Additional  assumptions  are  introduced  for  each  single  sub- 
domain,  and,  for  sake  of  clarity,  they  are  presented  in  the  specific 
sections. 

4.1.  Conservation  and  constitutive  equations 

Mass,  chemical  species,  momentum  and  charge  conservation 
laws  are  applied  in  the  different  domains.  Mass,  species  and 
momentum  conservation  equations  are  considered  in  the  gas  chan¬ 
nels.  The  charge  conservation  is  applied  in  the  anode,  the  cathode 
and  the  electrolyte. 

In  the  porous  media,  mass  transport  is  assumed  to  occur  mainly 
by  diffusion,  therefore,  momentum  balance  coupled  with  the  con¬ 
stitutive  law  for  porous  media  (e.g.  Darcy’s  law  or  Brinkman’s 
equation)  is  neglected.  This  assumption  has  been  verified  by  pre¬ 
liminary  simulations,  where  Brinkman’s  equation  is  implemented 
and  a  negligible  effect  is  demonstrated.  Figs.  3  and  4  present  the 
profiles  of  convective  velocity  field  in  a  section  of  anode  and  cath¬ 
ode  respectively,  when  Brinkman’s  equation  is  implemented.  As 
can  be  noted,  anode  and  cathode  convective  velocity  are  almost 
zero  (order  of  magnitude  10-5  m  s-1 ),  thus  convective  flows  (order 
of  magnitude  10-7  to  10-6  molnrr2  s-1)  in  the  electrodes  are  neg¬ 
ligible  compared  to  diffusive  flows  (order  of  magnitude  10-2  to 
1 0-1  mol  m-2  s-1 ).  The  diffusion  of  chemical  species  is  modelled  by 


Fick’s  law.  Convection  and  diffusion  phenomena  have  been  consid¬ 
ered  to  model  the  balance  of  chemical  species  in  the  gas  channels. 

The  charge  conservation  equation  is  applied  to  electrode  (elec¬ 
tric  charge)  and  electrolyte  domain  (ionic  charge).  The  constitutive 
law  to  describe  these  charge  transport  phenomena  is  Ohm’s  law. 

4.2.  Gas  channels 


According  to  the  conservative  and  consecutive  equations 
referred  in  Section  4.1,  the  set  of  equations  describing  fluid 
behaviour  in  gas  channels  is: 

+  V(pu)  =  0 

(3) 

+  u^(pYj)  =  V(pD, j-Vy,) 

(4) 

+  p(uVu)  =  -VP  +  uV2u 
ot 

(5) 

Eq.  (4)  is  obtained  from  the  conservation  law  of  chemical  species 
and  Fick’s  law,  considering  that  no  reactions  take  place  in  the  gas 
channels  (assumption  1).  Eq.  (4)  becomes  four  equations  when 
applied  to  the  chemical  species  in  gas  channel  streams:  H2  and 
H20  in  the  anode  side,  02  sand  N2  in  the  cathode  side.  In  obtaining 
Eq.  (5),  body  forces  within  the  gas  channels  have  been  neglected. 
The  assumption  of  ideal  gas  introduces  the  expression  of  density, 
which  completes  the  system  of  equations: 


The  dynamic  viscosity  of  a  gas  mixture  has  been  calculated  with 
Herning  and  Zipperer  formula,  proposed  by  Shan  et  al.  [18].  For 
the  single  chemical  species  of  the  mixture,  dynamic  viscosity  has 
been  taken  from  an  empirical  correlation  proposed  by  Todd  and 
Young  [1 9].  Fuller-Schettler-Gidding  correlation  [20]  has  been  used 
to  calculate  the  diffusion  coefficient  Dy  of  a  binary  mixture. 

4.3.  Electrodes 

In  this  case,  the  porous  nature  of  the  electrodes  is  taken  into 
account  in  defining  the  binary  diffusion  coefficients.  The  following 
expressions  are  derived,  respectively,  for  anode  and  cathode: 

=  V(pD02-N2,porousVyi)  (7) 
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Fig.  5.  Detail  of  the  single  fuel  cell  test  rig  in  the  fuel  cell  laboratory  at  the  University 
of  Perugia. 


dpYj 

dt 


V(/?Dh2  0-H2 ,  porous  ^  Yi ) 


(8) 


If  the  Knudsen  number,  Kn  =  A/Dpore,  is  near  1  or  higher  than  1, 
Knudsen  diffusion  is  important.  A  in  Knudsen  number  expression 
represents  the  mean  free  path  of  the  gas.  As  stated  in  [20,21  ],  Knud¬ 
sen  diffusion  should  be  considered  when  dealing  with  SOFC  porous 
electrodes  mass  transport.  The  diffusion  coefficients  in  the  porous 
media  take  into  account  Knudsen  diffusion  coefficient  DiK  which 
considers  the  interaction  between  molecules  and  pore  walls  and 
the  binary  diffusion  coefficient. 


D02,eff  = 


DN2,eff  = 
DH2,eff  = 


1  1 

- +  ■ 


%  porous  ^0  2,I< 


A 


1  1 

- + 


ij,  porous 


Au,/c 


A 


1  1 

- + 


ij,  porous 


Al2,/< 


DH20,eff  = 


1 


A 


ij,  porous 


Dh2o,k 


(9a) 

(9b) 

(10a) 

(10b) 


The  effective  diffusion  coefficients  of  a  binary  mixture  in  a 
porous  medium  are  calculated  through  the  following  expression, 
taken  from  Aguiar  et  al.  [22]: 

£ 

Aj,  porous  —  Aj  ~  (11) 

Besides  the  mass  transport,  the  electric  problem  is  described 
within  the  electrodes  domain.  The  charge  conservation  equation, 
coupled  with  Ohm’s  law  yields: 

^+V(aelV0)  =  O  (12) 

Additionally,  charge  conservation  equation  is  solved  to  deter¬ 
mine  the  current  density  distribution  within  the  electrodes  domain. 


4.4.  Electrolyte 


The  only  phenomenon  considered  in  this  domain  is  the  ion 
transportation,  O2-,  from  cathode  to  anode  interface. 


^Pion 

dt 


+  V(crV0ion)  —  0 


Table  2 

Model  parameters  and  material  properties. 


Parameter 

Units 

Calibrated  value 

Literature 

Fa, an 

[J  mol-1  ] 

128,000 

100,000l8l 

Fa, cat 

[J  mol-1  ] 

122,000 

120,000l8l 

Kan 

[Am-2] 

1  xlO9 

5.5  x  1 0s  181 

Kcat 

[Am-2] 

1  xlO9 

7x  lO8^ 

(7  Jon 

[^-1  m-1] 

2.26383 

2.26383^ 

Oan 

[ft-1m-1] 

30315.4 

30315.419! 

CTcat 

[ft-1  m-1] 

12792.4 

12792.419! 

S 

[-] 

0.3 

0.3™ 

T 

[-] 

10 

6llsl 

F>n2 

[m2  s-1  ] 

5.692  xlO-6 

8.322  xlO-6'18-21! 

Do2 

[m2  s-1] 

5.571  x  10-6 

8.322  xlO-6'18-21! 

Dh2 

[m2  s-1] 

2.416  x  10-5 

3.765  xlO-5'18-21! 

F)h2o 

[m2  s-1  ] 

1.407  x  10-5 

3.765  xlO-5'18-21! 

V</>ion  =  ^  (14) 

°ion 

5.  Boundary  conditions 

Unless  explicitly  indicated  otherwise,  internal-interface  bound¬ 
ary  conditions  are  set  to  ensure  continuity  between  two  contiguous 
domains,  x  and  y  directions  refer  to  the  directions  in  Fig.  1. 


5.1.  Gas  channels 

Boundary  conditions  for  Eqs.  (3)-(5)  for  both  gas  channels  are: 

•  No-slip  condition  of  the  gas  on  the  electrode/gas  channel  inter¬ 
face  and  gas  channel/current  collector  walls. 

•  Gas  velocity  entering  gas  channel  domain  is  known. 

•  Gas  channel  discharge  pressure  is  the  outlet  operating  pressure. 

•  Concentration  entering  the  gas  channel  for  every  component,  i, 
of  the  gas  mixture  is  known. 

•  For  channel  outlet,  only  convection  condition  is  considered. 

•  The  remaining  gas  channel  wall  is  isolated  for  mass  transport. 


5.2.  Electrodes 

Electrodes  are  in  contact,  on  one  side,  with  the  gas  and  the  cur¬ 
rent  collector,  and,  on  the  other  side,  with  the  electrolyte,  while  the 
remaining  surfaces  are  insulated. 

The  internal  boundary  conditions  to  solve  mass  transport  equa¬ 
tion  are: 


^  i 
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Fig.  6.  Simulated  and  experimental  polarization  curves. 
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Table  3 

Simulation  vs.  experimental  parameters. 


Voltage  [V] 

Uflab  [%] 

Ufsim[%] 

euf  [%] 

Jlab  [Am-2] 

jsim  [Am-2] 

ej  [%] 

0.998 

1.7 

1.6 

-5.88 

200 

190 

-5.00 

0.935 

8.7 

8.8 

1.14 

1000 

1025 

2.50 

0.885 

15.7 

15.8 

0.63 

1800 

1839 

2.16 

0.839 

22.7 

22.8 

0.44 

2600 

2650 

1.92 

0.795 

29.7 

29.7 

0 

3400 

3449 

1.44 

0.751 

36.7 

36.5 

-0.54 

4200 

4244 

1.04 

0.705 

45.5 

43.3 

-4.83 

5000 

5036 

0.72 

0.656 

50.7 

49.8 

-1.77 

5800 

5788 

-0.20 

0.596 

57.7 

55.9 

3.11 

6600 

6469 

-1.60 

0.498 

64.7 

64.1 

-0.92 

7400 

7444 

0.59 

0.456 

66.5 

67.5 

1.2 

7600 

7821 

2.90 

Fig.  7.  Result  for  the  current  density  in  a  2D  domain  of  a  disk-shape  anode  supported 
SOFC(V=  0.656  V). 

•  at  the  electrode/electrolyte  interface,  gas  species  gener¬ 
ated/consumed  by  electrochemical  reactions  (1)  and  (2),  are 
related  to  the  current  density  by  the  following  equations: 

SJVh2  =  (15) 


MOLAR  FRACTION  -  ANODE  GAS  CHANNEL 


Fig.  8.  Molar  fraction  of  fuel  components  in  the  interface  anode/gas  channel  along 
the  cell  (V=  0.656  V). 


Fig.  9.  Hydrogen  and  water  mole  fractions  variation  inside  the  anode  (V=  0.656  V). 


^H20  =  2f  ^h2o 

(16) 

n-(jd) 

n'No2  =  4F  M02 

(17) 

The  boundary  conditions  to  solve  the  electrical  problem  within 
electrodes  are: 

•  at  electrodes/current  collector  interfaces,  the  boundary  condition 
is  given  by  the  electrical  potential 

|  JeUn  =0  (lg) 

0el,cat  —  ^cell 

Conditions  ( 1 8 )  are  applied  under  the  assumption  that  the  poten¬ 
tial  drop  within  the  current  collectors  is  negligible,  which  is 
justified  by  the  high  electrical  conductivity  of  the  current  col¬ 
lectors,  compared  to  those  of  the  electrodes. 

•  at  electrode/electrolyte  interfaces,  the  boundary  condition  for 
the  electric  problem  is  given  by  the  generated  electrical  current 
defined  with  the  Butler-Volmer  equation: 
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nJ=Jo 


(  nfiFr)  \ 

eXP(^-J-eXP 


-n(l  - P)Fruct^ 


RgT 


(19) 


where  Jo  is  the  exchange  current  density  and  is  computed  for  anode 
and  cathode: 


,  _  y  PH2  (PH2o\ 

J0,an-yan—^  —  ) 


-0.9 


exp 


— £act,a 


RgT 


Jo, cat  —  /cat 


mc 

PrefJ 


exp 


-£act,a 

RgT 


(20) 


(21) 


r] act  represents  the  activation  loss,  and  is  the  difference  between 
the  ideal  and  real  potential  at  the  electrode/electrolyte  interface: 


?7act  —  ^electrode  ^electrolyte  ^rev 


(22) 


1/rev  is  computed  with  Nernst’s  equation,  for  anode  and  cathode 
[11]: 


Vr£ 


F  _j_  in  (  PH20 
|  =  _£°  +  _ln/ 


(23) 


CURRENT  DENSITY 


RgT  . 

In 

nF 


(24) 


Note  that  expressions  (23)  and  (24)  assume  that  most  of  the 
potential  contribution  is  provided  by  the  anode.  This  assumption 
is  justified  by  the  fact  that  the  electrochemical  reaction  takes  place 
at  the  anode  side. 


5.3.  Electrolyte 

There  is  no  generation  of  ionic  current  inside  the  electrolyte 
domain  but,  as  stated  in  assumption  1,  in  Section  4,  at  the  elec¬ 
trolyte/electrode  interface.  Thus,  an  ionic  current  flows  through 
the  boundaries  in  contact  with  anode  and  cathode  while  the  other 
boundaries  of  the  electrolyte  are  insulated. 

From  the  stoichiometry  of  reactions  (1)  and  (2),  two  electrons 
pass  through  the  external  circuit  when  generating  one  oxygen  ion. 
Thus,  ionic  current  is  a  half  of  electronic  current. 


Fig.  10.  Current  density  distribution  for  several  voltages. 


6.  Model  validation  and  results 

The  set  of  equations  defined  in  Sections  4  and  5  are  implemented 
and  solved  for  steady-state  conditions,  using  the  commer¬ 
cial  software  for  Finite  Element  Method  resolution,  COMSOL 
MULTIPHYSICS®. 

Due  to  the  symmetry  of  the  cell,  it  is  assumed  that  each  radial 
section  behaves  in  the  same  manner,  thus  a  2D  representation 
of  this  section  is  used  as  the  cell  domain  (Fig.  1).  Operational 
parameters  used  for  the  simulation  are  those  used  during  a  set  of 
experimental  tests  conducted  at  the  Fuel  Cell  Group  of  the  Univer¬ 
sity  of  Perugia  (Fig.  5).  Operating  conditions  are:  1073  K,  ambient 
pressure,  humidified  hydrogen  for  the  anode  feed  (97%  H2  and  3% 
H20)  and  synthetic  air  as  the  oxidant  flow  (21%  02  and  79%  N2  from 
industrial  purity  gas).  Inlet  gas  flows  are  24  N1  h-1  for  the  anode  and 
60  N1  h- 1  for  the  cathode. 


xio4 
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Fig.  12.  Current  density  distribution  vs.  molar  fraction  of  H2,  in  anode  side,  and  02, 
in  cathode  side. 


In  order  to  obtain  the  anodic  and  cathodic  flow  velocity  in  the 
considered  2D  section,  the  following  relation  is  applied: 


Qan,c 


2^rmanifhChannei 


Qan.cat  — 


(25) 

(26) 


It  should  be  noted  that  the  choice  of  solving  the  problem  in  a 
2D,  rather  than  a  3D  domain,  leads  to  neglect  the  channel  width 
variations  along  the  radius,  and,  as  the  final  consequence,  to  ignore 
the  velocity  reduction  of  the  gas,  as  this  proceeds  along  the  gas 
channel.  The  present  study  is  limited  to  obtain  preliminary  results 
on  a  2D  domain,  and  to  verify  the  electrochemical  performance 
predictions,  and  the  convergence  speed  of  the  entire  numerical 
simulation  process.  It  must  be  stressed  that  the  present  preliminary 
results  provide  relevant  useful  indications  on  the  cell  operational 
conditions,  with  an  acceptable  accuracy  of  the  results. 

When  comparing  the  calculated  ideal  open  circuit  voltage  (OCV) 
to  the  experimental  OCV,  the  introduction  of  an  empirical  factor, 
ofocv.  is  required: 


c^ocv  = 


OCVmeasured 

OCVcajcuiatecj 


(27) 


Several  reasons  explain  the  discrepancy  between  the  calculated 
and  measured  OCV.  Specifically,  for  the  disk-shape  SOFC,  Momma 
et  al.  [23]  measured  a  back-diffusion  of  cathodic  gas  into  the  outer 
region  of  the  anode.  However,  the  investigation  of  such  phenomena 
is  beyond  the  scope  of  the  present  paper. 

Experimental  works  conducted  at  Forschungszentrum  Jiilich, 
Germany  [24]  and  unpublished  results  obtained  at  Postech  Univer¬ 
sity,  in  South  Korea  [25],  show  that  when  assembled  and  operated 
in  an  entire  fuel  cell,  single  components  show  different  electrical 
proprieties,  compared  to  those  of  single  samples  not  embedded 
in  a  fuel  cell.  Both  [24,25]  noticed  a  conductivity  reduction,  not 
directly  correlated  to  the  anode  and  cathode  materials  themselves, 
but  to  the  contact  between  them  and  the  electrolyte.  As  a  conse¬ 
quence  of  these  effects,  data  available  in  the  literature  for  SOFC 
components  conductivity  are  not  directly  applicable  for  fuel  cell 
simulations.  For  this  reason,  these  values  have  been  extrapolated 
from  the  comparison  with  the  experimental  data. 

By  comparing  the  results  of  the  model  with  the  experimental 
data,  a  calibration  of  the  parameters  is  obtained.  The  calibrated 
parameters  and  the  corresponding  values  found  in  the  literature 


3  ACTIVATION  LOSSES  -  ANODE 
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Fig.  13.  (a)  Activation  losses  at  anode  interface  and  (b)  activation  losses  at  cathode 
interface. 


are  compared  in  Table  2.  Fig.  6  illustrates  that  the  experimental  and 
calculated  polarization  curves  are  mostly  coincident  but  in  the  ini¬ 
tial  part,  where  losses  related  to  energy  of  activation  are  dominant. 
Table  3  presents  the  values  of  voltage  versus  current  density  and 
coefficient  of  utilization  of  the  fuel  (simulated  and  experimental) 
and  the  corresponding  percentual  error.  The  error  never  exceeds 
5%  for  the  current  density,  and  reaches  5.8%  for  the  fuel  utilization. 

Fig.  7  shows  the  distribution  of  the  current  density  in  the  entire 
domain  when  the  cell  voltage  is  0.656  V.  The  current  density,  as 
expected,  diminishes  while  hydrogen  is  consumed  along  the  cell. 

Fig.  8  presents  hydrogen  and  water  mole  fractions  profiles  along 
the  cell  length  in  the  gas  channel  (interface  gas  channel/electrode), 
when  the  cell  voltage  is  0.656  V.  Fig.  9  shows  H2  and  H20  mole  frac¬ 
tion  variation  within  the  whole  anode  domain,  and,  together  with 
Fig.  8  illustrate  the  impact  of  the  occurrence  of  hydrogen  oxidation. 

Results  shown  in  Figs.  7-9  are  referred  to  a  specific  electric  load 
demand,  characterized  by  a  cell  voltage  of  0.656  V.  This  value  can 
be  considered  as  a  typical  cell  voltage  for  SOFC. 
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IONIC  POTENTIAL  ANODE-  ELECTROLYTE  INTERFACE 
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Fig.  14.  (a)  Reversible  voltage  at  anode-electrolyte  interface,  (b)  electronic  potential  distribution  at  anode-electrolyte  interface,  and  (c)  ionic  potential  distribution  at 
anode-electrolyte  interface. 


However,  it  is  important  to  investigate  the  internal  condition  of 
a  single  cell  for  other  load  requirements. 

As  illustrated  in  Fig.  1,  since  the  current  collection  is  performed 
via  a  metallic  mesh  located  on  the  entire  surface  of  the  electrode 
in  the  gas  channel,  the  current  path  is  perpendicular  to  the  cell 
plane,  i.e.  the  current  density  is  uniform  along  the  cell  thickness, 
thus  one  dimensional  graph  can  fairly  describe  the  current  density 
distribution  inside  the  cell. 

The  current  density  profile,  for  several  applied  voltages,  is 
shown  in  Fig.  10.  As  the  applied  voltage  decreases,  the  total  cur¬ 
rent  provided  by  the  cell  increases.  Since  the  inlet  fuel  and  oxidant 
flow  rates  are  constant,  the  fuel  and  oxidant  utilizations  increase 
as  the  applied  voltage  decreases.  For  high  fuel  or  oxidant  utiliza¬ 
tions,  a  gas  starvation  takes  place  locally,  in  some  regions  of  the 
two  electrodes,  thus  generating  a  strong  current  density  reduction. 
As  visible  in  Fig.  10,  for  a  voltage  as  low  as  0.556  V,  the  current 
density  reduction  along  the  cell  is  very  pronounced. 


As  shown  in  Fig.  11,  for  voltages  below  0.5  V,  the  outer  part  of 
the  cell  is  not  properly  exploited,  i.e.  the  current  density  becomes 
close  to  zero  (75  times  smaller  than  maximum  local  value  of  cur¬ 
rent  density)  at  the  exit  of  the  cell.  Eqs.  (1 5)-(l  7)  show  that  current 
density  is  directly  related  to  hydrogen,  water  and  oxygen  con¬ 
sumption/production,  thus  gas  starvation  and  the  drop  of  current 
density  are  strictly  connected.  Such  operational  condition  should 
be  avoided,  especially  for  long  term  operation,  because  it  leads  to 
severe  materials  degradation  issues  [26-28]. 

As  shown  in  Fig.  12,  in  the  present  case,  the  current  density 
reduction  is  due  to  oxygen  starvation  at  the  cathode. 

A  better  distribution  of  the  current  density  and,  consequently,  a 
reduction  of  local  gas  starvation  could  be  achieved  via  an  increasing 
porosity  of  the  electrodes  from  the  inlet  to  the  outlet.  With  the 
present  fuel  cell  technology,  however,  it  is  preferable  to  operate  the 
fuel  cell  at  reduced  fuel  and  oxygen  utilization.  It  should  be  noted 
that  when  a  stack  is  operated,  the  stack  temperature  is  regulated 
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3  REVERSIBLE  POTENTIAL  CATHODE-  ELECTROLYTE  INTERFACE 


IONIC  POTENTIAL  CATHODE-  ELECTROLYTE  INTERFACE 


Fig.  15.  (a)  Reversible  voltage  at  cathode-electrolyte  interface,  (b)  electronic  potential  distribution  at  cathode-electrolyte  interface  and  (c)  ionic  potential  distribution  at 


cathode-electrolyte  interface. 

by  the  air  flow  rate.  The  energy  balance  imposes  an  air  flow  rate 
that  is  generally  within  a  range  of  3-7  times  the  stoichiometric 
value  [29,30].  For  this  reason,  oxygen  local  starvation  is  not  an  issue 
for  SOFC  when  operated  in  a  stack.  On  the  other  hand,  in  order  to 
achieve  high  conversion  efficiency,  fuel  utilization  needs  to  be  kept 
as  high  as  possible.  An  engineering  solution  is  to  recycle  the  anode 
outlet  to  the  inlet,  so  that  the  overall  fuel  utilization  of  the  system  is 
kept  high,  even  if  the  fuel  cell  operates  with  reduced  fuel  utilization 
[31-33]. 

Fig.  13a  and  b  illustrates  the  activation  losses  along  the  cell 
length.  According  to  assumption  1,  Section  4,  the  electrochemical 
reactions  take  place  at  the  electrode/electrolyte  interface;  there¬ 
fore  activation  losses  occur  in  these  surfaces  (reduced  to  lines  in 
the  2D  schematization). 

Activation  losses  depend  on  gas  concentration  and  current 
density,  as  shown  in  relations  (19)-(24).  The  high  oxidant 
utilization  corresponding  to  V=  0.556  V  leads  to  an  extremely 
low  oxygen  concentration  and,  as  a  consequence,  to  differ¬ 
ent  distributions  of  activation  losses  for  the  anode  and  for  the 
cathode. 


Lenght  [m] 

Fig.  16.  Trend  of  ionic  potentials  at  anode  and  cathode  interfaces  for  low  voltages. 
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OHMIC  LOSSES  -  CATHODE 


Fig.  17.  (a)  Ohmic  losses  in  anode  domain,  (b)  Ohmic  losses  in  cathode  domain,  and  (c)  Ohmic  losses  in  electrolyte  domain. 


Fig.  14a-c  presents  the  variation  of  the  terms  used  for  the  calcu¬ 
lation  of  the  anode  activation  losses  in  order  to  identify  the  cause 
of  this  behaviour.  Considering  the  orders  of  magnitude,  it  may  be 
stated  that  the  anode  activation  losses  are  controlled  by  the  dif¬ 
ference  between  reversible  potential  and  ionic  potential  at  the 
electrolyte-anode  interface  (electronic  potential  is  negligible),  Eqs. 
(22)  and  (23). 

Fig.  15a-c  shows  the  terms  involved  in  calculation  of  the  cath¬ 
ode  activation  losses,  Eqs.  (22)  and  (24).  In  this  case,  the  electronic 
potential  is  not  negligible  and  its  evolution  along  the  length  of 
the  cell  is  mostly  constant,  thus  the  distribution  of  the  cathodic 
activation  losses  is  mainly  controlled  by  the  ionic  and  reversible 
potentials. 

The  effect  of  oxygen  starvation  for  V=  0.556  V  or  lower  volt¬ 
ages  is  directly  connected  to  the  difference  of  ionic  potential 
between  anode-electrolyte  and  cathode-electrolyte  interfaces. 
Fig.  16  shows  the  trend  of  these  two  potentials  at  both  interfaces 
for  low  voltages  (0.556  V  and  0.456  V)  due  to  progressive  starvation 


of  oxygen  near  the  outlet.  In  this  region,  it  is  clearly  visible  that  the 
difference  of  potentials  for  0.456  V  becomes  zero,  leading  to  zero 
current  density. 

If  a  PEN  lumped  structure  is  assumed  for  modelling  the  fuel  cell 
[4-9],  the  mass  transport  inside  the  porous  media  is  not  directly 
modelled.  Flowever,  the  gas  concentration  in  the  bulk  is  different 
from  that  at  the  reaction  zone,  due  to  mass  transport  phenomena, 
thus  the  voltage  calculated  with  gas  concentration  in  the  bulk  is 
higher  than  that  calculated  at  the  reaction  zone  (TPB).  As  a  con¬ 
sequence,  numerical  models  assuming  a  PEN  structure,  usually 
introduce  so-called  concentration  losses,  which  quantify  the  volt¬ 
age  reduction,  due  to  the  mass  transport  phenomena  within  the 
electrodes.  This  loss  is  deduced,  together  with  the  activation  and 
the  ohmic  losses,  to  OCV,  in  order  to  calculate  the  cell  voltage  [4-9]. 

Although  with  the  modelling  approach  adopted  for  the  present 
study,  it  is  not  necessary  to  introduce  such  a  computation  for  eval¬ 
uating  the  cell  voltage,  it  is  useful  to  evaluate,  as  a  post  processing 
operation,  the  effect  of  the  mass  transport. 
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This  is  done  by  calculating  the  difference  between  Nernst  volt¬ 
age  at  the  bulk  and  at  the  reaction  zone: 


RgT . 

?7conc,an  —  — ppr  In 


P' 

1 


Ho 


^H20,bulk 


?7conc,cat  —  — 


RgT , 

^|n 

nF 


H2,bulk 


02 

D02,bulk 


H20 

1/2 


(28) 


(29) 


is  operated  in  an  isothermal  atmosphere,  controlled  by  an  elec¬ 
tric  furnace.  An  extension  of  the  model  will  be  to  implement 
energy  equations  and  to  compare  the  resulting  temperature 
distribution  to  that  of  a  single  cell  embedded  in  a  stack. 

(3)  The  implementation  of  the  model  at  non-steady-state- 
conditions,  and  in  a  3D  domain. 
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As  expected,  since  the  anode  is  much  thicker  than  the  cathode, 
the  anode  concentration  loss  is  much  higher  than  the  cathodic  one. 
This  trend,  however,  changes  for  very  high  current  densities  and 
fuel  utilizations,  i.e.  for  a  low  cell  voltage. 

Finally  Fig.  17a-c  represents  the  ohmic  losses  at  the  anode, 
cathode  and  electrolyte,  respectively.  Since  the  electrical  and  ionic 
conductivity  are  assumed  homogenous  within  each  domain,  the 
trend  in  ohmic  losses  is  the  same  of  the  current  density. 

7.  Conclusions 

Numerical  simulations  of  SOFC  are  still  in  an  early  stage  of  devel¬ 
opment  because  of  the  relative  novel  interest  in  SOFC  and  the 
complexity  of  the  process  coupling. 

Most  of  the  current  mathematical  models  are  not  capable  to 
simulate  the  entire  domain,  because  the  assumed  lumped  PEN 
structure.  In  addition,  they  are  restricted  to  just  one  specific  geom¬ 
etry. 

The  detailed  model  presented  in  this  paper  is  based  on  partial 
differential  equations  and  can  be  applied  to  different  geome¬ 
tries  and  cell  materials.  In  the  present  study,  the  model  was 
implemented  and  validated  on  a  disk-shape  planar  single  cell. 
Experimental  data  for  validation  are  obtained  at  the  FC  lab  of  the 
University  of  Perugia.  Numerical  and  experimental  results  show  a 
good  agreement,  with  an  error  below  5%  for  the  current  density, 
and  6%  for  the  fuel  utilization.  The  results  obtained  allow  visualiza¬ 
tion  of  the  evolution  of  the  most  important  physical  variables  in  the 
entire  fuel  cell  domain,  thus  identifying  limiting  phenomena,  such 
as  local  gas  starvation  at  high  fuel  and  oxidant  utilization,  and  the 
contribution  of  mass  transport,  ohmic  losses  and  activation  losses 
to  the  voltage  reduction. 

8.  Further  development 

Three  main  limitations  of  the  model  are  amenable  to  further 
developments  to  create  a  more  robust  and  efficient  tool  for  design¬ 
ing  SOFC  and  predicting  the  behaviour. 

(1)  This  model  refers  to  humidified  hydrogen  as  fuel.  To  simulate 
the  behaviour  of  a  fuel  cell  operating  on  other  gases,  a  math¬ 
ematical  model  capable  of  describing  internal  reforming  and 
equilibrium  (water-shift)  reactions  within  the  anode  side  must 
be  defined. 

(2)  In  this  paper,  energy  balance  is  not  taken  into  consideration 
because  the  fuel  cell  from  which  experimental  data  are  taken 
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